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Abstract. In this paper we discuss three symbolic approaches for the generation of a finite 
difference scheme of a partial differential equation (PDE). We prove, that for a linear PDE with 
constant coefficients these three approaches are equivalent and discuss the applicability of them 
to nonlinear PDE's as well as to the case of variable coefficients. Moreover, we systematically 
use another symbolic technique, namely the cylindrical algebraic decomposition, in order to 
derive the conditions on the von Neumann stability of a difference scheme for a linear PDE with 
constant coefficients. For stable schemes we demonstrate algorithmic and symbolic approach 
to handle both continuous and discrete dispersion. We present an implementation of tools 
for generation of schemes, which rely on Grobner basis, in the system Singular and present 
numerous examples, computed with our implementation. In the stability analysis, we use the 
system Mathematica for cylindrical algebraic decomposition. 
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Introduction 



The finite difference method for hnear PDE's befongs to the very classical topics in mathematics. 
However, its exposition in the classical books like [24] often contains mysterious steps, relying 
on the huge experience, gathered in last centuries. An algebraist is often confused with such 
exposition and asks, whether there is a way to split the whole picture into purely analytic and 
algebraic parts and how is it possible to automatize the process of scheme generation and further 
analysis of its properties. Hence the ideas to generate finite difference scheme in an algebraic (or 
a symbolic) way are folklore, see for instance [H [19] for approaches and older implementations. 
Terminologically, we address a difference scheme as symbolic polynomial expression involving 
unknown function and shift (or difference) operators. We do not attach initial and/or boundary 
conditions to differential resp. difference equations, since the generation of a difference scheme 
as above is independent on them. As we will see, von Neumann stability can be seen as global 
result, always being a necessary condition for stability of a problem with initial and/or boundary 
conditions (and sometimes sufficient condition as well). Of course, one uses initial and/or boundary 
conditions for numerical solving, but the splitting of the whole problem into purely symbolic pre- 
processing and numerical post-processing seems to be the way to address such problems in the 
future. 

In the article [11 , Gerdt et al. used for the first time several new ideas like the use of integral 
relations in discrete form (especially useful if one deals with conservation laws), the formulation 
of the scheme generation problem as a task for difference elimination and the systematic use of 
involutive and Grobner bases. Inspired by these ideas, we present our approaches, which will make 
the overall picture of scheme generation and analysis more complete. 

The ideas of algebraic analysis suggest the separation of a problem into analytic and algebraic 
parts. This allows, in the case of linear PDE's, to treat systems of equations via modules over 
algebras (D-module theory, homological algebra etc.). In such a case there exist many algorithms 
and several powerful implementations. Grobner and involutive bases play a fundamental role in 
such algorithms, see e. g. [22 1 120 1 13]. 

In the other direction, the differential algebra (e. g. [H]) and difference algebra (e. g- [5]) theories 
allow one to tackle nonlinear equations as well, though the algorithms in these realms are very 
complicated. In particular, up to now we do not know any implementation of a basis construction 
algorithm for the difference algebra. Notably, a new algorithmic approach to (nonlinear) difference 
equations seem to follow from the letterplace approach [IT]. However, one needs to elaborate all 
details of this promising direction. 

The usage of the famous cylindrical algebraic decomposition (CAD) originated from real alge- 
braic geometry to von Neumann stability problems goes back to [191 [S] . Since that times more 
implementations of the CAD evloved and their performance has been greatly enhanced. 
This paper is organized as follows. We start with minimal prerequisites and revisit the basic 
concepts of scheme generation, paying attention to the algebraic background including Grobner 
bases and elimination tools in the Section [I] 

We discuss three symbolic methods, used in applications in the Section [2] and prove their equiv- 
alence in the case of a linear PDE with constant coefficients in the Theorem 12.21 In cases of 
other PDE's only one method will work in general, see Remark 12.31 As for linear PDE's as in 
the Theorem, we propose to use module formulation and Grobner bases for eliminating module 
components, which can be seen as a natural generalization of the Gaussian elimination to matrices 
over rings. We show the merits of this method, applied to the classical equations of mathematical 
physics (heat, wave, advection equations) for various approximations. Note, that the method we 
propose can handle high order approximations, which are seldom used in the theory of PDE, but 
quite often in the theory of ODE as high order Runge-Kutta methods. 

In the Section [4] we present an algebraic and constructive formulation of von Neumann stability 
via ring homomorphism. We shortly revisit the concepts of cylindrical algebraic decomposition 
and connect its use to the questions, arising from difference schemes. 
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In the Section [5] we consider the A- wave equation Utt = Xu^^ and perform both generation and 
stabihty analysis of several difference schemes, obtained with different approximations. We demon- 
strate the merits of the semi-factorized form of a difference scheme, it turns out to be especially 
useful for higher dimensional situation. 

In the Section [6] we demonstrate, that the determination of continuous respectively discrete dis- 
persions for a PDE respectively its difference scheme can be algebraized to the large extent as 
well. 

All the examples in this paper have been computed with our implementation of tools for difference 
scheme in a freely available computer algebra system Singular jT5]. The corresponding library 
f indif s . lib will be distributed with the next version of SINGULAR. For the cylindrical algebraic 
decomposition we use the commercial system Mathematica; indeed there are freely available 
systems like QEPCAD and REDLOG, which are able to do the decomposition as well. 
From the viewpoint of applied mathematics there is a general skepticism about the use of symbolic 
methods. With this paper we want to stimulate a discussion between scientists of both fields based 
on a realistic viewpoint. 

1. Algebraization of differential and difference equations 

1.1. Types of operator algebras. First of all we have to fix a computable field k, our base 
field, it is mostly the field of rational numbers Q or complex rational numbers Q[i\. (Computing 
with reals or complex numbers is in principle possible, but only with a fixed precision, i.e. with 
a rational approximation, or one can compute with algebraic extensions in roots of polynomials 
- this being not so interesting for our purpose.) We can extend the base field by indeterminate 
constants, i.e. rational functions in the constants: K = k{a, b,c, . . .), which can be specialized to 
special values at any step of our computation if necessary. K is called the field of constants. 
We fix a set of variables x := (xi, . . . , a;„) and an 'algebra' of functions C = C{x) in the x's, for 
instance differentiable functions or functions in discrete (shiftable) arguments. C is not our object 
of computation. Instead we consider various operator algebras, consisting of operators, which act 
on C. There are many operators, which one can handle in this framework, for example 

a) multiplication with a variable: Xi : u{x) € C i— >■ Xiu(x) G C; 

b) multiplication with a function f ^ C: nif : u{x) G C i— > f{x) ■ u{x) G C. 

c) partial differentiation: di : u{x) G C h^ q^^' G C; 

d) partial shift operators: Ti : u{x) G C h^ u(xi, . . . ,Xi + I, . . . , x„) G C; 

e) partial A-shift operators: T^ : u{x) G C H> u{xi, . . . ,Xi + X, . . . , x„) G C, clearly Ti ~ T}; 

f) partial difference operators: Ai = T/' — 1; u{x) G C M- T^{u{x)) — u{x) G C; 

g) g-dilation operators: Dd : u{x) G C i-^ u{xi, . . . , qxi, . . . , Xn) G C, 
h) et cetera. 

Fix a set S of operators, we consider the operator algebra A :— K{S) being the sub algebra of 
all (linear) operators HomK{C,C) generated by S, i.e. the smallest linear subspace closed under 
multiplication of operators. As long as S consists of a finite number of pairwise commuting and 
independent operators the resulting algebra is isomorphic to a polynomial ring: _ftr[ii, . . . , t,„]. 
Otherwise we get (non- commutative) quotient algebra of the free algebra K{S) by the two-sided 
ideal of all relations of S. 

Example 1.1 (Algebras with constant coefficients). The algebras of linear partial differ- 
ential and shift (or difference) operators with constant coefficients are commutative K-algehras, 
isomorphic to K[xi^ . . . , x„] . We denote them by K[di, . . . , 9„] and K\Ti^ . . . , r„] respectively. 

Example 1.2 (Algebra with polynomial coefficients). The algebra of linear partial differ- 
ential operators with polynomial coefficients is the Weyl algebra. It is non-commutative but has 
simple commuting relations. We denote this algebra as K{x,d \ dx — xd + 1) , what means that 
this algebra is generated by the {x, d} over the field K , and has the multiplication, defined on its 
generators: x ■ d = xd, d ■ x = xd -\- 1. The generalization to multivariate case is easy, the variable 
{xi} commutes with the variable {dj} expect for the case j = i, then the relation as above applies. 
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Why the variables satisfy this relation? Consider the multiplication of two operators, x := Xi 
and d := ^-. Take some differentiable function f and apply the Leibnitz rule to the product: 
{dx){f) = d{xf) = xd{f) + f = (xd + !)(/)• Hence, in the operator form {dx — xd ~ 1)(/) = 0. 
Consider the algebra of linear \- shift operators with polynomial coefficients, having in mind 
A — Ax. Like before, we ca derive the relation between operators T :— T^ and x. For any function 
in discrete argument f, {Tx){f) = T{xf) = T{x)T{f) = {x + Ax)T{f) = ixT){f) + Ax ■ T{f). 
Hence, in the operator form this relation becomes Tx ~ xT + Ax ■ T. The algebra, corresponding 
to the difference operator A = T — 1 has the relation Ax — xA + A + 1. These algebras are 
so called G-algebras, in which Grobner basis algorithms exist and are implemented in the system 
Singular:Plural (\T3\). see e. g. [I8] . 

Example 1.3 (Algebra with coefRcients in rational functions). Algorithmic computations 
are possible in the algebras whose coefficient fields are rational functions in x: K{x){d \ dx — 
xd + 1) and K{Ax,x){T \ Tx — xT + Ax ■ T), which are called rational Weyl algebra resp. 
rational shift algebra. Algebraically .speaking, a passage from polynomial algebra to a rational 
algebra may be achieved by means of localization. 

Example 1.4 ( Differential and Difference Algebra). In order to handle non-linear differen- 
tial resp. difference equations with polynomial nonlinearities, one can consider a full differential 
resp. difference algebra K[{0^u \ /3 € N"}], where O stands for differential resp. difference oper- 
ators. Note, that O^u is a variable, representing 0^{u), where u = u{x) is a dependent variable 
and X an independent one. Note, that these algebras are commutative and their infinite generating 
sets are algebraically independent (but dependent differentially resp. in a difference way). 
The given nonlinear equations can be taken as generators of the differential resp. difference ideal 
(that is an ideal, closed under the action of corresponding operators) in an above algebra. 
Since such algebras are infinitely generated, they are not Noetherian. Grobner basis-like algorithms 
are therefore not terminating in general. Nevertheless some parts of the theory from the linear 
situation is extended to this general situation. 

In this paper we work algorithmically with linear partial differential operators with constant 
coefRcients. However, in some parts we address and discuss more general situations as well. 

1.2. Presentation of a system of differential equations. Any partial differential equation 
(depending of its kind) defines an element of an algebra of corresponding differential operators. A 
solution of a system of equations fulfill any equation of the left ideal generated by the equations of 
the system in the algebra. Hence, the solution does not depend on the choice of a basis (that is, a 
generating set) of the ideal. The first possibility of applying symbolic algebra is to compute a better 
basis of the ideal, like a Grobner basis or an involutive basis (like Janet basis, see |20[ [Q] 123) - as 
far as it is possible. The advantage of thus a pre-processing could be: check the consistency of the 
system of equations, find hidden constraints or integrability conditions of the system, determine 
the dimension of the solution space etc. 

These data are well known for standard equations from mathematical physics, but the methods we 
propose are methodologically applicable to any system of equations. Let us recall a small example 
(by W. Seller [23]) as an illustration. 



Example 1.5. 



Mz + y Wx == ( Uy,, + yuxy + u:c =0 ^ „ 

Uy ^0 [ Uxy = Wy2 = 



Hence, the initial system is equivalent to {ux = Uy ^ Uz =0}. 
In the case of a linear system (S) 

n 

Si = y Dij • Uj, i — 1, . . .m, Dij G A 

we associate to (S) the submodulc P = P{S) C A" generated by the columns of a presentation 
matrix D G Mat{m,n;A), and, finally, a factor-module M{S) := A"/P{S). We can simplify 
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the system finding a special presentation matrix, or we can read properties of the system firom 
computable invariants of the module M{S). 

In the example above, the system can be written as I ^ q ) • ("") = ( n ) ■ Hence, the system 

algebra is A = K{x, y, z){dx, dy, dz \ dxX = xdx + 1, dyy ~ ydy + 1, dzZ — zdz + 1) (that is the 3rd 
Weyl algebra) and the presentation matrix for the system module M(S), written in columns of the 
original presentation (that is, transposed to the usual raw presentation) is P{S) = {dz+ydx, dy) £ 
A^^^. As a submodule of A, P{S) is an ideal and it has two polynomial generators {dz + ydx, dy}. 
The Grobner basis of P(5) is equal to Q{S) = {dx,dy, dz) and hence, M{S) = A/Q(S) ^ K{x, y, z) 
as A-module. Thus, d[mK(x,y,z) M{S) = 1 is the dimension of the solution space of S. 

1.3. Grobner basis algorithm and elimination tools. The notion of Grobner basis can be 
given in a common way for different classes of algebras. Recall the basic notation for monomials 
and monomial ordering. We shall use the short notation 9" := 9"^ 9^^ . . . 9"", a S N". Finitely 
generated operator algebras, which we are dealing with, have infinite dimension as ii'-vector spaces. 
The infinite set of monomials constitutes this basis. 

• For operator algebras in operators {9i, . . . ,9„} with constant coefficients, the monomials 
are {9" | a G N"}, they form the basis of an algebra over the field K. 

• In the case where the coefficients are polynomial in {xi, . . . , x^}, the monomials are 
{x" ■ d^ \ a E N™, /3 G N"} and they form the basis of an algebra over K. 

• When the coefficients are rational functions, the monomials {9" | a G N"} constitute the 
basis of an algebra over K{xi, . . . , x„i). 

We dealing not only with ideals of an algebra A, but also with submodules of a free module 
A"^ — (Bi^iAci, where e^ is the canonical i-th basis vector. We extend the notion of a monomial 
to A^ by supplying a monomial with one of the basis vectors. Clearly, if Mon{A) := {nia} is the 
set of monomials of A, bijective to N*", then the monomials of A^ are {niaei \ a G N", 1 < i < r}. 

Definition 1.6. A (global) monomial ordering on an algebra A as before is a total ordering -< on 
the set of monomials Mon(A) bounded from below and compatible with the multiplication, i.e. it 
fulfills the the following conditions for all a,/3,7 G N".' 

• 1 ^ rua, 

• ma -< mfj => niamy -< mpm^. 

Since -< is total, any nonzero polynomial f (z A can be uniquely sorted according to its monomials. 
The highest term (that is monomial times nonzero coefficient) is called the leading monomial of 
f. We say, that rUa \ rn/j (ma divides rnjs), if\fl<i<nai<l3i. Note, that the divisibility is a 
partial ordering. 

Any monomial ordering is extendable to a module monomial ordering in several ways. The most 
common ways are: either sorting module monomials first by the monomial ordering and then by 
the number of component, or first by component and then by the monomial ordering. 

Definition 1.7. Given a monomial ordering -< of A, then monomial orderings <top (term-over- 
position) and ^pot (position- over-term) on the set of monomial of A'' are defined by: 

{ma, ei) <top {mp,ej) iff {ma -< mp or if ma = mp, then i < j), 

respectively 

{ma,ei) <pot {mi3,ej) iff {i < j or if i = j, then ma -< m^s). 

Definition 1.8. A Grobner basis of a submodule M d A^ is a finite subset G d M that satisfies 
the following property. For any / G M \ {0}, there exists a element of the basis g G G, such that 
the leading monomial lm{f) = uiaCi is divisible by lm{g) ~ mpCi, i. e. j3 \ a. 

An immediate but very important application of Grobner basis is the normal form of a vector of 
polynomials. Namely, if G = {gi, . . . ,gm} is a Grobner basis of a submodule M C A^ , then for 



6 V. LEVANDOVSKYY AND B. MARTIN 

any v e A'' there exist w, a^ e A, such that 

m 

V — y aiQi + w, where either atgi = or hn(w) ^ \Ta{aigi) and eitherw = or lm(it;) ^ lm(w). 

4 = 1 

One denotes w = NF(u, G) and calles uj to be a normal form of v with respect to G. Note, that 
v€ M = {G) ii and only if w; = NF(w, G) = 0. 

A Grobner basis G = {gi, . . . , gm\ is called reduced^ if for any 1 < i < rn and j 7^ i, no monomial 
of 5j is divided by lm((7i). Any Grobner basis can be made reduced in a finite number of steps. 
Notably, normalized reduced Grobner basis (that is, having I's as leading coefficients of gis) is 
unique. It turns out, that normalized reduced normal form is unique, hence we address it by 
saying the normal form. 

There are effective ways to compute a Grobner basis, like the Buchberger's Algorithm, involutive 
algorithm and Faugere 's F^ or F5 algorithm. Grobner bases have been implemented in all major 
computer algebra systems. More details for the commutative case can be found meanwhile in 
any standard textbook on computer algebra, e. g. in [M] . Consult with the 21 HH] for the non- 
commutative case of operators with variable coefRcients. 

It is important to mention, that applied to a module, generated by the columns of a constant 
matrix, the result of a Grobner basis algorithm (with respect to position-over-term ordering) is 
identical to the result of Gaussian elimination. 

Using the special monomial ordering, a Grobner basis algorithm can be used to eliminate some 
of the variables {ui | z e /} of a given system, i.e. to compute a basis of Mj := Af n Aj, 

A^ Ai{ui,ie I). 

Lemma 1.9. (Elimination of variables). Let -< be an elimination monomial ordering for 
{ui \ i (z 1} on Mon(j4) (that is, ma G Aj and j ^ I implies ma -< Uj). Let G be a Grobner basis 
of M , then G f] Ai is a Grobner basis of Mj. 

Note, that a lexicographical ordering of monomials by ui > U2 > ■ . • > Wn induces an elimination 

monomial ordering for any set {ui, . . . , w„}. 

We can also eliminate module components, and usually it is much easier than the elimination of 

variables. 

Lemma 1.10. (Elimination of components). Let G be a Grobner basis of a submodule M C 
A^ with respect to the module monomial ordering -<pot, let Fg := Aei © • • • © Acs C A^ be the free 
submodule of the first s components, then G H Fs is a Grobner basis of M H Fg. 

The proofs of both lemmata on elimination are easy and can be found in e. g. [141 118) for various 
situations. 

Remark 1.11. We want to stress the fact, that the algorithm we propose, using the opera- 
tor formulation of the problem, will always be faster, than the algorithm in difference algebra, 
which is used by Gerdt et al. in [llj . The difference in complexity lies in the number of vari- 
ables/components one uses and the intrinsic differences between two similar-looking elimination 
concepts. Computing in the case, when the functions and discretizations of their derivatives 
u,Ut,Uxx etc. appear as variables (difference algebra approach), one must distinguish between 
the two multiplications, namely the action one (denoted by •), appearing when difference opera- 
tors act on u 's, and the operator one (denoted by ■), used for composition of difference operators 
themself. One either uses the involutive basis approach with its division of all variables into mul- 
tiplicative and non-multiplicative, or forbids the multiplications between u 's in the Grobner basis 
algorithm. In addition, one has to employ a complicated elimination ordering, which respects the 
special role of u 's. 

By passing to the submodule over a ring of difference operators, we do not use u 's at all. We use 
the linearity of operators in order to present the equations as linear operators themselves, presented 
by polynomials in difference operators with constant coefficients, involving parameters. Thus, we 
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use less variables, simple and efficient module component elimination orderings. Moreover, we 
shift the attention of Grobner basis algorithm from single polynomials rather to their components, 
which results in easier and faster computation, not speaking on much optimized usage of memory. 

2. Three Equivalent Approaches and the Main Theorem 

Assume we are dealing with m spatial variables xi, . . . , x„ and one temporal variable t = a;„i+i. 
We denote x" := x"^ ■ ■ ■ x^t"'"+^ for a S N^+i and \a\ — ^ a^. Then we use notations 

dx°' [[ "^i 

A single linear PDE with constant coefBcients can be written as follows 

(1) ^ - X! '^f^^^^ = ^' 

fieB 

where B C N"'+^ is a finite set and cp ^ K. 

In order to compute approximations, we rewrite the equation in an arbitrary inferior point t — 

(ii, . . . , im, n) e G := ZA.ti x • • • x ZAx„ x ZAi (where G is often identified with Z™+-^) of the 
grid, that is 

(2) Y. ^^<^ = 0- 

On the grid, one needs to give an approximation to a function u'' f, through a linear combination of 
similar expressions, coming from differentials of order, lower than /3, that is Vt € G, V/3 (z B \ {0} 

(3) Afj = u^^fs — 2_] d^jul-, with | 7 |<| /3 \,d-yj € K. 

7erj6G 

Note, that all but finite number of d^j are zero. 

We say, that a general problem of approximation of a partial differential equation is well-defined, 

if for any /3 G _B \ {0} (that is for any u^/i , appearing in the equation with non-zero coefficient 
except for u itself) there is a unique approximation from A/^. 

Indeed, on the grid G we have natural shift operators Tj,. : u' M- w'"'''^', where e^ is the i-th 
canonical basis vector. That is T^,. (w(. . . , i^, . . .)) = w(. . . , t^ + Ax^, . . .)) is a well-known forward 
shift operator. Clearly T^^ is invertible operator, since the inverse is just a backward shift operator. 
Thus, working with T^- we allow exponents to be integers. For an exponent vector a G Z™"*"^, 
denote T" = T^"/ • . . . • T^"'"r""+\ In what follows we wih use the field of fractions K{T) := 
K{T,„. . .,T^STt) of the polynomial ring K[T] := K[T^,, . . . ,T,^,Tt]. 

Thus, it is possible to use shift operators and rewrite the previous equation in a single generic 
point t of the grid, as the following lemma shows. 

Lemma 2.1. In the notations as before, there exist exponent vectors S{l),6{j), J G G such that 

(4) T^'-'^u^^y = J2 d^jT^^^Ku^-yf or {u^of = ^ d^jT^^^^-^^'\u.^.Y . 

7erjeG 7erjeG 

Proof. For 1 < fc < m -I- 1 set k^ :— m\n{]k,i'k | .7 G G,7 G ^^d^.j 7^ 0}. Then, setting 
5{t) := I — K and respectively S{j) := j — R we obtain two exponent vectors for the monomials in 
shift operators. D 

According to the lemma, in what follows we will derive and encode approximations for functions 
on the grid with the help of shift operators. This allows us to drop the grid point notation as soon 
as shift operators are present. In other words, we use shift operators to formulate the problem in 
a generic point of the grid. 
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There are several approaches to the computation of a finite difference scheme of a single partial 
differential equation. 

2.1. Mimicking Difference Algebra Approach. Consider the formal consequences of equali- 
ties P as in dl]) and A/3 as in ^ over the commutative ring Rb ■= K{T)[u^i3 \ (3 € B]. Recall, 
that the variables {u^ii} are algebraically independent. In other words, we consider an ideal / of 
the ring Rb, generated by P U {Ap | /3 S P}. Since B is finite, the ring Rb is Noetherian and 
contains a subring R := K{T)[u]. Hence the ideal J :— I Ci K{T)[u] exists and it is computable 
(e. g. by the elimination of all but one variables as in Lemma Tl.Qp . Since P is a principal ideal 
domain, J is generated by a single element, say p € K(T)[u\. Clearing denominators, we obtain 
a polynomial expression p £ K[T][u\. Dividing its leading coefficient out, we obtain the unique 
result. 

Note, that we do not work with difference ideal in a difference ring, but just mimic a technique 
for presenting the objects. 

2.2. Algebraic Analysis Approach. Let us order the set {u^/i \ (3 £ P} according to the 
monomial ordering (hence by the total degree as well), starting, say, from the highest appearing 
Pmax, and write the resulting ordered list as a column vector U. Then it has the form U = 
[u^fima.^ , ■ ■ ■ , u]'^ ■ Since P and Aj^ are linear equations with coefficients in K[T] in the entries of U, 
we put each equation as a row in a matrix, say M, with entries in iir[T], such that MuU — 0, where 
• stands for the action of shift operators with coefficients in K on functions in discrete arguments. 
Then we can perform operations from the left on the matrix M only, without engaging the u's. 
This can be seen as the application of one of the principles of algebraic analysis. We compute then 
the intersection of iir[T]-module M with the free submodule, generated by the last component of 
U. The latter is an ideal J C K[T] of all polynomials p in shift operators, such that pu u — 0. 

2.3. Term Rewriting System Approach. Consider the equations from (j4|) in the monic form, 
that is 

Let us treat them as rewriting rules for symbols {u^(> \ /3 £ B \ {0}}, which rewrites every 
appearance of W3./3 with the sum on the right hand side. We call such system 5*. Since the latter 
always involves Ux-t if and only if 7 ^ /3, we do the following. At first, we order occurring variables 
as in the approach l2.2| that is descending with respect to the monomial ordering {u^iimax , ■ • ■ , u\. 
Then, in the same sequence, the variable u^'^ is substituted (replaced) with the right hand side 
of the corresponding approximation A^. According to the ordering, the result of the substitution 
does not contain variables, which are higher, according to the monomial ordering. Hence the result 
of the substitution is unique and consists of shift operators, acting on u only. 

Theorem 2.2. For a single linear partial differential equation with constant coefficients, the 
following three methods lead to the same result and hence they are equivalent: 

a) As in \2.1[ I D K{T)[u\ = (/), where f is chosen to be monic polynomial from K[T][u\; 

b) As mm for r = \B\ holds K[TY D Mn K[T]er ^ J C K[T]. Moreover, J = (g) , where 
g £ K[T] is chosen to be monic; 

c) As in \2.S\. P -^s h and h is normalized in addition. 

Proof, a) Ap is already a Grobner basis in K{T)[ux(>] by the product criterion, because the leading 
monomials of its elements are coprime (since the set {u^/i \ (3 G B} is algebraically independent. 
Since NF(P, Ap) e K{T)[u], we obtain that {NF(P, Ap)} U A^ is a Grobner basis of P U Ap. The 
uniqueness follows from the uniqueness of normalized reduced normal form [14. 

b) Proceeding with the vector as above and starting with higher leading monomials, the matrix 
representation M of the set Ap is already upper triangular (with its entries in /'ir[r]). Moreover, the 
last row has at exactly two nonzero elements (say, the last two ones in that row). Thus M is already 
in a row-reduced form. Note, that making complete reduction of the rows will produce a matrix 
M', where each row contains exactly two nonzero elements: (0, . . . , 0, /i(P), 0, . . . , 0, fr{T)) with 
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fi, fr G I'^[T]- Hence A/' simplifies the set of approximations and correspond to the completely 
reduced Grobner basis. 

Now, we append to M (or M' , what is equivalent) the row P', corresponding to the discretized 
equation P. The computation of the Gaussian elimination of the resulting matrix amounts in 
reductions of P' with the rows of row- reduced M'. The result of such reductions is a row vector 
with the only nonzero entry f{T) G K[T] at it last position. Since it is the same as the reduced 
normal form NF(P', M') and M' is a Grobner basis, the ideal J of the Theorem is principal ideal, 
generated by /(T). Then the normalized /(T) is single operator, acting on u, that is the difference 
scheme. 

c) Note, that the application of rewriting rules in the sequence, as described above, leads to the 
NF(P, A^). Since by a) Ap is a Grobner basis with respect to any monomial well-ordering, this 
normal form is unique. Fixing a monomial ordering, we can produce another rewriting system S' 
by applying rewriting rules to every right hand side of S, starting ascendingly from the smallest 
nonzero /3 e -B. Then S' becomes {{u^^ -)■ fp{T)u) \ /? e S\{0}}, where fp{T) e K[T^^] C K{T). 
It is straightforward, that S" does not depend on the sequence of reductions like S anymore and 
S' is confluent. The reduction of P with respect to S" is the same as with respect to S, hence it 
is unique. D 

Remark 2.3. Note, that the equivalences of the previous Theorem, do not hold in general. Namely, 
Algebraic Analysis Approach \2.2\ works only for linear PDE, since it relies on the module struc- 
ture, which is linear per definition. Moreover, both a) and b) do not necessarily work in the case 
of variable coefficients. In the latter, there are different concepts for discretization. As soon as 
one deals with algebras, where x and T^ do not commute anymore, the left normal form of a 
vector (the computation uses subtractions of left multiples of an approximation and thus invokes 
non- commutative multiplication) is not necessarily the result of rewriting of any term u^fi (which 
just plugs the right hand side expression into the place where the term resides, just not invoking 
non-commutative multiplication) . 

Indeed the only method, which works for all cases of systems, is the Term Rewriting System 
Avvroach \2.3\ (Ap is a Grobner basis for arbitrary linear approximations) . 

3. Generation of difference schemes 

Armed with the methods from the previous section, we proceed with the generation of schemes 
for linear equation with constant coefficients. We prefer the method of algebraic analysis from l2.2l 
in contrast to Gerdt et al. [TT] , who used the method of difference algebra because of the reasons 
of practical complexity, which is significantly lower in the approach 12.21 However, Gerdt et al. 
systematically follow the difference algebra approach for nonlinear equations, where in |11) they 
have obtained an interesting nice-behaving scheme of the order 3 for the the original equation of 
order 2, which however does not contain switches as usual schemes. 

A big class of PDE's from mathematical physics are of conservative type in space-time with 
coordinates {t,x) G M"*: 

dP/dt - dQ/dx = 0. 
By Gauss formulae this equation can be transformed into its conservation law form: 



/ {Qdt + Pdx) = 0. 

JdV 



IdV 

We choose some discretized integration contours and approximations rules for the integrals and 
proceed as above. The difference schemes, which we obtain by elimination are by construction 
fully consistent [TT] . 

3.1. Approximation Rules and their Operator Form. The most general way to approximate 
a PDE is to use the integral relations like J "^^ utt{x,t)dt — ut{x,tn+i) — ut{x,tn) together with 
further approximations of differentials and integrals. 
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However, a large class of equations might be written in a so-called conservation law form, which 
can be obtained e.g. by applying the Green's formula. For example, the equation -^ ^ ^ = is 

equivalent to the equation f Pdx + Qdy — for arbitrary piecewise smooth closed contour T . 

r 
Contour approximations. There are many possibilities for choosing contours and approxima- 
tions. Since we are using rather rectangular than quadratic grids, the two most used approximation 
on contours are the node points of the rectangle and the middle points of a grid with a double 
distance, as illustrated by the pictures below. 



+ i 


r " 


1 1 


I. I 






- 1 


■ 1 


1 ► 



J-l 



J + 1 



fc-1* 



J-1 



.7 + 1 



Although by applying the Green's formula we lower the order of an equation by 1, the approxima- 
tion formulas, derived from the contour are usually more complicated, than the approximations, 
derived from the original equation and integral relations. Note, that this is not a problem for an 
implementation, since complicated manipulations with polynomial expressions can be performed 
effectively with modern computer algebra systems. 

Approximation of differentials via Taylor series. Applying the Taylor expansion up to the 
2nd order, we obtain u{x ± A.t) — u{x) ± AxUx{x) H — §-Uxx{x) + 0{Ax^). Hence, 
u{x + Aa;) — u(x) 



Ux{x) 



Aa 



0{/\x) (forward difference) 



U(x) — U(X — Ax) ^^/ A \ /1 1 1 1 rr- 

or Ux[x) = h 0(Ax) (backward dmerence). 

Ax 

Subtracting these two equalities we get u{x + Ax) — u{x) + u{x) — u{x — Ax) = 2AxUx{x) +0 (Ax^) , 

hence Ux{x) — ^^ 2a" — ~^ 0{Ax'^) (central 1"* order difference). 

Adding these two equalities and rewriting the result, we get the central 2"'' order difference 



u{x + Ax) — 2u{x) + u{x — Ax) 
Ax^ 



Uxx{x)+0{Ax^). 



Approximation of integrals. Closed Newton-Cotes formulas give rise to e.g. trapezoid an 
pyramid rules, whereas open Newton-Cotes formulas lead us to e.g. midpoint rule. 



The trapezoid rule is expressed as follows: 

xq-\-Ax 



f{x)dx - iAx(/(a;o) + f{xQ + Ax)) - -^f" {£), x^ < i < x^ + Ax. 



Hence, in the approximation we obtain, for f{x) — Ux{x): 

X{) + /\X 



u{xq + Ax) - u{xq) = / Ux{x)dx = -Ax{ux{xn) +Ux{xn + Ax)), 



and hence {Tx — 1) • u — ^Ax{Tx + 1) • Ux. 
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Pyramid (or Simpson's) rule looks as follows: 

xa+2Ax 

f{x)dx = \^x{f{xo) + 4/(a;o + As) + /(xq + 2Ax)) - ^^/(4)(^), 

xa 

hence its difference form is | Aa; • (T^ + 4Tx + \) • u^ — {T^ — 1) • u. 
Open Newton-Cotes formula for 1 point, 

a;o+2Aa:; 

fix)dx = 2Ax/(a;o + Ax) + 0{Ax^f'), 

xo 

leads us to the midpoint formula Ax ■ T^ • Ux = {T^ — 1) • m. 

Summary. In a form of difference operators, we gather the most used approximations 

• Forward diflference {Ax, 1 — T^) • {ux,u)'^ = 

• Backward difference (Ax ■ T^, 1 ~ T^) • {ux, u)^ — 

• A 1"* order central appr. (2Aa: • T^, 1 - T'^) • {ux,uf = 

• A 2"'* order central appr. (-Ax^ • T^, (1 - Txf) • {uxx,uf = 

• Trapezoid rule {\Ax ■ (T^ + 1), 1 - Tx) • {ux,uY = 

• Midpoint rule (2Ax -Tx,!- T^) • (ux, uf = 0. 

• Pyramid rule (i Aa; • (T^ + 4^^ + 1), 1 - X?) • {ux, u)^ = 

• Lax methocfl {2At ■ Tx, T^ - 2TtTx + 1) • (ut, u)^ = 

• Parametric temporal difference for < < I: (At- (6^* + (1 -6*)), l-Tt) • (ut,u)^ = 0. 
li 6 — resp. = 1, it becomes forward resp. backward difference. 

Assume that the difference scheme involves quantities Axi, . . . , Axm, At and originates from a 
typical set of approximations. By definition, difference scheme is of the smallest difference order, 
hence the shift polynomial p, describing it, is irreducible. It turns out, that in many situations we 
want to present p as the sum of products of operators. Taking into account the future applications 
to the von Neumann stability, we propose the following. 

Definition 3.1. A semi-factorized presentation of a linear difference scheme of order 0{Axi , . . . , Ax\^ , Af^ 
to he the sum p = Ax^pi + . . . + Ax^Pm + Af^pt for pi e K[T], such that pi does not involve 
Axi in its coefficients and most (if not all) pj do not involve Axi, . . . , Xm, At. 

Unlike nodal form, a semi-factorized form allows compact descriptions of very complicated and 
higher dimensional schemes. Note, that in the examples we present, it turns out that there exists 
a unique (up to constant factors) semi-factorized presentation. In the implementation we have a 
method for computing a semi-factorized form constructively. 

Example 3.2. Consider the ID heat equation ut — c?Uxx ~ for a parameter a. 

We approximate ut with backwards difference At ■ Tt • Ut = (Ti — 1) • u, or, in the nodes of the 
grtd,At-{u,)T+'^{u)T+'-{u)T. 

As for Uxx, it is approximated with the 2nd order weighted centered space method, that is 

Ax^ ■ Tx • Uxx = {OTt + (1 - 0)) ■ {Tx - 1)^ • u, where < 9 < I. 
In such a way, we obtain a matrix formulation of the problem 

1 -a' \ /"* 

-At-Tt Tf-l • I Uxx I = 0. 

-Ax^TxTt {9Tt + {1 - 6)) ■ {Tx - l)y 



used in e.g. discretization of the advection equation 
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By computing a Grobner basis (with the algebraic analysis approach), we obtain a single polynomial 
in shift operators for the scheme 

-a'^AteT^Tt+a^At{9-l)T^+{2a^Ate+Ax^)T^Tt-{2a^At{e-l)+Ax'^)T^-a^AteTt+a^At{e-l) 

Its semi-factorized form is Ax'^Tx{Tt — l) — a'^At{Tx — l)'^{0Tt + l — 6) = 0. In the following example 
we show Singular code for obtaining these objects and for producing a nodal presentation of the 
scheme, which is 

^ • {u]tl u-^,) ^ ■ {u-Xl 2U-XI + «r^) - ^(""+2 - 2<+i + uj) =. 0. 

The scheme we obtained is called FTCS ifd = 0, BTCS if 9 = 1 and Crank- Nicholson, if 6 = ^. 

If is easy to see, that this scheme is consistent with the original differential equation for any G M. 
Namely, since "^+^'^"^+^ = w* + 0(At) and "^+' ~^2^'''^"^"^' = u,^x + 0{Ax'^), we have 

—ut - eu^^ - (1 - G)u.^.^ = —ut - u.^x == ©(At) + 0{Ax^). 
a'' a^ 

Thus, the order of the scheme is (Ai, Ax^). 

Example 3.3. In this example we do computations with SINGULAR and findifs. lib. As we see 

from the matrix formulation, our parameters are At, Ax, a, 0, to which we add the parameter d, 

which will be needed later for the check of stability. The variables of our ring are Tt and T^ . We 

define the ring in SINGULAR and the matrix of equations as follows: 

ring r = (0,a,dx,dt ,theta,d) , (Tx,Tt) , (c,Dp) ; 

matrix M[3] [3] = 

1, -a~2, 0, // the equation itself 

-dt*Tt, 0, Tt-1, // appr. u_t with backward difference 

0, -dx~2*Tt*Tx, (theta*Tt+(l-theta) ) *(Tx-l) ~2; // appr. u_xx with theta-centered space 

Now we transpose the module and then call the std routine for getting the Grobner basis. 

module R = module (transpose (M) ) ; module S = std(R) ; 

print (S) ; 

==> 0, 0, 1, 

0, (-a~2*dt)*Tt, (-a~2) , 

S[3,l],Tt-l, 

As we can see, the first column vector is the only one with the values in the 3rd component only. 
S[3, 1] is displayed since this polynomial (which describes the difference scheme) is big. 
poly p = S[3,l]; p; // assign and print the answer 

==>(-a"2*dt*theta)*Tx~2*Tt+(a"2*dt*theta-a~2*dt)*Tx~2+(2*a'2*dt*theta+dx'2)*Tx*Tt+ 
(-2*a~2*dt*theta+2*a'2*dt-dx~2)*Tx+(-a~2*dt*theta)*Tt+(a'2*dt*theta-a'2*dt) 

We proceed with the semi-factorized form and visualization. 
LIB "f indif s . lib" ; // load the library for schemes 
ideal I = decoef (p,dt) ; // auxiliary procedure 
I; // the sum of elements of I gives p 
==>! [1] =(dx'2) *Tx*Tt+(-dx'2) *Tx 

I[2]=(-a~2*dt*theta)*Tx~2*Tt+(a'2*dt*theta-a~2*dt)*Tx~2+(2*a'2*dt*theta)*Tx*Tt+ 
(-2*a~2*dt*theta+2*a'2*dt)*Tx+(-a~2*dt*theta)*Tt+(a'2*dt*theta-a'2*dt) 
From this structure, we can obtain the semi-factorized operator form of the scheme: 
f actorizeCl [1] ) ; // we suppress the output 
f actorizeCl [2] ) ; // factors with multiplicities 
==>[1]: 

_[l] = (-a-2*dt) 

_[2]=Tx-l 

_ [3] = (theta) *Tt+(-theta+l) 
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[2]: 

1,2,1 

Hence, the semi-factorized form is Ax'^T^iTt — 1) — o^AtiT^ — l)^{9Tt + 1 — 9) = 0. 

list L; L[l] = theta; 

dif poly2tex(I ,L) ; // we show though only a part of this string 

==>\frac{-l}{a--[2} \tri t}\cdot (u-{n+l}_{j+l}-u-{n}_-[j+l>)+ ... 

The string above (in tex format) is the nodal presentation of the scheme, which has already been 
demonstrated in the previous example. 

4. Symbolic methods for von Neumann stability analysis 

4.1. Stability rings, morphisms and polynomials. We refer the reader to e. g. [24] for details 
about stability. Suppose that t is the temporal variable and xi, . . . , Xm are the spatial variables. 
We start with a finite difference scheme, written in the nodal form on the uniform orthogonal grid 
with steps At, Axi, . . . , Axm- We suppose to work in the interior region, which is bounded, say, 
by Li, . . . , Ljn in spatial directions. 

In the von Neumann stability analysis, one presents the functions on the grid as discrete Fourier 
modes, that is 

m 

fe=i 

where % is a linear map, g is a new symbolic variable, < IkAxk < L^- We abbreviate 
f^jk '■— T^lkAxk- Then, one substitutes this presentation of nodes into the equation, performs 
simplifications and obtains a polynomial G in one variable g with constant coefficients. 

The von Neumann stability criterion (as in [24]) states, that if for every root ^ of G, one has 

1^1 < 1, then the difference scheme is stable. 

The Lax equivalence theorem can be stated in the following form (adopted from [53]). A 
consistent scheme for a well-posed linear initial value problem is convergent if and only if it is sta- 
ble. For a well-posed linear initial-boundary-value problem, however, stability is only a necessary 
condition for convergence in general. 

We do not address algorithms for algorithmic check of consistency of a difference scheme with its 
differential equation in this paper. There are several methods for doing it using algebraic tools 
like [8] [12] . However, for several equations we treat we show the usage of semi-factorizcd form of 
a scheme to the positive conclusion about such consistency. 

Let A be the algebra of functions on a given grid. It is naturally a module over the algebra R 
of linear partial difference operators with constant coefficients G[Tt,Tx-^ , . . . , Tx^] over some field 
G 3 Q(Ai, Acci, . . . , Axm)- The action of R on discrete Fourier nodes, using the map Xi can be 
written as follows: 

X(^ • <n..,J - 9" ■ X{ul,,..,j and x(lt « ^.n..,J = ^^'^'^^^^'^ ' x{ul,,..,J for all j,. 
The map x and this action give rise to the constructive homomorphism of C = Q(At, Axi, . . . , Aa;^)- 
algebras 

X ■■ G[Tt,T^^,. . .,T^J — > G{[i,sinx^,coSx^,. . . ,sin^^,coSj;^]/J„i)[g], 

where J^ — {i^ + 1, sin^.^ + cos-^^ — 1, ... , sin-^^^ + cos'^^^ — 1) is the ideal. We denote this con- 
structive stability morphism by the same letter x and note its C-linearity. It is defined by its 
values on the generators of the source algebra and is given by x{Tt) — 9 and xi^j^) — e''='^^^= — 
cos (3s + i ■ sin /3s for all 1 < s < to. 

The constructiveness of this approach and hence, its applicability in computer algebra systems, 
lies in the following. We choose the basic numeric field to be exact complex-rational numbers 
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Q[i]/ {i" + 1) . On demand, we can do further algebraic extensions of this field. Moreover, we avoid 
complex exponentials by passing to the sine and cosine, incorporating their natural algebraic re- 
lations in the factor ideal. Then, a stability morphism can be implemented and used in every 
computer algebra system, being able to compute Grobner bases. 

Now, let P — J2a a '^a,aTiT" be the operator form of the finite difference scheme P»u = 0, where 
T^ stands for T;^} ■ T^^ for a multi-index a = (ai, . . . , a™) G N". Then, 

m 

X(P) = Y.^-^o.x{Ttrx{T,r - E^a.o HicosPk + t -sm /Skr^g'^ = E^-5" 

a^a a,Oi fc— 1 a 

is the univariate polynomial in g, which we call the stability polynomial of a given difference 
scheme. Obviously, the degree of x(-P) is the same as the highest degree of Tt in P. 



Example 4.1. Let us continue with the examvle \3.2l In order to prepare the scheme for stability 
analysis, one can rewrite it as 



y»+i „," 



i' - «;Vi = ^'d{0 ■ (u-+^ 2u]+l + u]+') + {i-e)- (u;v2 - 2u;vi + «")), 



with d :— — — -■ We prefer to work with the semi-factorized operator form of the scheme Ax'^TxiTt— 

1) - a^At{Tx - ifiOTt + 1 - 0) = 

Creating the stability ring and performing simplification and factorization in it (see next example 

for the Singular code), we obtain the following linear polynomial in the variable g. 

(i cos -I- sin) • {{{-2a^de) sin +2a^de + I) ■ g + {2a^de ~ 2a^d) sin ~2a^de + 2a^d - 1) 
The first factor i ■ cos(/3) + sin(/3) = e^'^ is ignored in stability analysis, since it is of magnitude 1. 

Example 4.2. We continue with the Examvle \3.3[ Define the semi-factorized scheme again. 

poly P = Tx*(Tt-l) + (-a'2)*d*(Tx-l)"2*((theta)*Tt+(-theta+l)) ; 
ring r2 = (0,a,theta,d) , (Tx,Tt) , (c,Dp) ; 
poly P = imap(r,P); 

Now, we create the stability ring ST (which will be Q{a, 9, d)[g, i, sin, cos]) and a map x from r2 
(it isQ{a,9,d)[Tx,Tt] to it: 

ring ST = (0,a,d,theta) , (g, i ,sin,cos) ,lp; 

ideal Rels = std(ideal(i2+l , sin~2+cos~2-l) ) ; // the ideal of relations 

map chi = r2, ideal (sin+i*cos ,g) ; 

poly P = chi(P); // the mapping 

P = NF(P,Rels); P; // reduction modulo the relations in Rel 

==>(-2*a"2*d*theta)*g*i*sin*cos+(2*a~2*d*theta+l)*g*i*cos+ . . . 
ideal FP = f actorize(P) ; // factorization 

The polynomial together with its factorization have been presented in the previous example. 

So, we came from a system of linear equations to a single univariate polynomial in the stability 
ring. The next problem we face is the following: 

Given a univariate parametric polynomial P, find out, under which conditions on 
parameters all the roots of P lie in the complex unit circle. 

As it have been already mentioned in [HI [1^ , this problem can be solved algorithmically with the 
help of the CAD (Cylindrical Algebraic Decomposition). 
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4.2. Cylindrical Algebraic Decomposition. The algorithm for CAD goes back to G. Collins 
et al. It is one of the most important algorithms, used for quantifier elimination not only in real 
algebraic geometry [1]. On the other hand, its algorithmic complexity is high and can be double 
exponential in the number of variables. Nevertheless, the universality of the method makes it very 
powerful and applicable to various problems. 

A finite set of polynomials {pi, . . . ,Pm} S M[xi, . . . , a;„] induces a decomposition (partition) of 
M" into maximal sign-invariant cells. A cell in the algebraic decomposition of {pi, . . . ,Pm} G 
M[a;i, . . . , Xn] is a maximal connected subset of M" on which all the pi are sign invariant. 

Definition 4.3. For n g N, let Wn : M" — > M"^^, (xi, . . . ,x„_i, a;„) H* {xi, . . . ,Xn-i) de- 
note the canonical projection. Let {pi,...,p„i} G (Q)[a;i, . . . , a;„]. The algebraic decomposition 
of {pi, . . . ,Pm} is called cylindrical, if 

• For any two cells C,D of the decomposition, the images 7r(C),7r(Z?) are either identical or 
disjoint. 

• The algebraic decomposition of {pi, . . . ,Pm} H Q[a;i, . . . , a;„_i] is cylindrical. 

For instance, any algebraic decomposition of M^ is cylindrical. 

There are several sophisticated implementations of the CAD algorithm. We are using the one from 
the system Mathematica, where two commands, CylindricalDecomposition and Reduce are 
available in the context of CAD. There are also freely available systems QEPCAD by C. Brown 
d] and REDLOG by A. Dolzmann et al. ^. 

4.3. CAD and Stability. 

Example 4.4. Let us continue with the examples \3.2\ [^TT] Let us represent the root of a stability 
polynomial as jr , where c = 2a^d(l — 9) sin— 2a^(i(l — 0) + \, d' — 2a^d9{l — sin) + 1 
Since d' > 0, we have to solve the inequality —d < c < d, that is c + d > and d > c. The first 
inequality 2a^d{l — sin) > is always satisfied, and the second is equivalent to a^d{20 — 1)(1 — 
sin) + 1 > 0. In this example, we compare the functions CylindricalDecomposition and Reduce 
of Mathematica 

CylindricalDecomposition [{a~2*d*(2*theta-l)*(l-s) + 1 >= 0, 
-1 <= s <= 1, a > 0, d > 0}, {theta, a, d, s}] 

returns 

(e < i&&a>0&&(0<d< ^ — )&&-l<s<l) I 

(d> -^ — && 1 - a d + 2a de ^^^ -^N 11 (0> i fc&a>o&&d>0&& -Ks < 1). 

^ -2o?+Aa?e -a?d + 2a?de - - > ^^ ^ 2 - - ) 

Whereas executing more specialized 

Reduce [a > && d > && <= theta <= 1 && 

ForAll[s, -1 <= s <= 1, a'2*d*(2*theta - l)*(l-s) + 1 >= 0] , {theta, d}] 

gives us more informative and structured answer 

a>Okk(0<e<-kkO<d< ^ — ) 11 (i < 6K 1 && d > 0), 

from which we conclude, that 

• if h 1^ S <1, the scheme is unconditionally stable 

• if Q < 9 < — , the scheme is stable under the condition d = -j-^ < — — -r. 

■' - 2 A^ - 2a^{l-29) 

The quantity d — -^-^ if often called Courant (or Courant-Friedrichs-Lewy) number. It is classical 
to express conditions on the von Neumann stability in terms of the Courant number. 

Example 4.5. Consider the ID advection equation ut -\-aUx = 0. We approximate ut with the 
parametric temporal method and Ux with the trapezoid rule. As a result, we obtain the difference 
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scheme in the semi-factorized form Ax ■ (T^ + 1) • {Tt — 1) + 2aAt • (T^ — 1) • (9Tt — {9 — 1)) = 0, 
which reads as follows in the nodal form: 

2^ • i-]ti «"+! + -r' <) + -^ ■ ((^i-iti -r) -(0- i)«+i - <)) - 0. 

^45 one can easily see, this scheme is consistent with its differential equation. The stability polyno- 
mial is linear with complex coefficients, so we present it as a fraction. The reformulated stability 
problem, which we have to solve, is 

4a^d^i29-l) ^TT 

-2 < ^ ,;. ... < 0, V/3 ^ -Z 

Since t := i-^-mi — 0? the right hand side inequality is equivalent to 9 < i The left hand side 
is equivalent to Aa'^d^{29 — 1) + 2{Aa'^d'^{9 — 1)^ + t)) > 0. Since t G [0,oo), we have to show 
that < 4a2d2(26l - 1) + 802^2(6* - l)^ = 402^2(^2 + (51 _ i)2)^ ^j^^t is true for all d. Of course, 
computations with CAD confirm this answer. 
Thus, this scheme is unconditionally stable if 9 < -^ o,nd unstable otherwise. 

5. Examples for A- wave equation 

We consider a parametric equation utt — ^'^Uxx ~ with parameter A 7^ and its higher dimensional 
versions. We construct finite difference schemes for several different approximations and analyze 
them for stability. 

5.1. Conservative law with parametric time approximation. The presentation via the con- 
servation law is f X^Uxdt + utdx — 0. We use trapezoid rule for both for the contour integral and 

r 
spatial integral relations. For temporal integral relations, we employ the parametric difference 

with 9 e [0, 1]. We obtain the following system of difference equations: 

^Ah-i-T,Tt + T,+Tt-l) X'^At-iT^Tt-Tt-T^ + l) \ fuA 

iAx • (T, + 1) 1 - tJ . U, = 

At ■ {9Tt + {1 - 9)) l-Ttj \uj 
After the computation of Grobner basis, we obtain the scheme 

° = 2Z^ • ("^"^' - 2"^"^' + "^") - 2^ • (<+2 - ^^^ + "^+2) + 
+ X^ • (^ • Kll - 2«;>r + <+') - (2^ - 1) • «:2 - 2U-+1 + u]+') + {9 -I)- («7+2 - 2^7+1 + u])) 

The stability polynomial of 2nd degree is rather complicated. However, factorization reveals a 
factor g — 1. The other factor is linear, but with complicated coefficients. We present it as 
g — ■^. Since both c and d' are complex numbers, we compute the absolute value of them. Then, 

lld'll = {4:9^d^ - 1) ■ {cos{(3^) - 1) - 2 and ||c|| = ||d'|| - 4d'*(20 - l)(cos(/3^) - 1). 

Then, ll^ll < 1 ^ < 4^^-(/3^/2)^p,^r3^^^^7^yT^ < l- 
Consider the left hand side inequality 

< 4^^ Bin(/3^/2)- (,,.,. j;;,:^;;^^,).^, - o < (2. - m^^d^ - DM^^r + 1) 

Since A9^d'^ > <^ A9'^d'^ — 1 > — 1 > -, tt^j the second factor is always positive. Hence, 

sin(/3a;/2)2 

the inequality is satisfied as soon as > -j- 



(2^-1) 



2- 

46*2^4 sin(;3^/2)2 + 1 - sin(;3^/2)2 > (4^^ sin(/3^/2)2)(26i - I) ^ 



The second inequality reads as Ad sn\(Bx/2) , „^ ,, ^ 7-; — 7—7; < 1. Then, 

^ -^ ^' ' ' (46*2^4- l)sin(/3^/2)2 + l - 
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^"^^^/'^' >(2.-i)^(^-if+.r.^t/t,>o, 



4d4sin(/3^/2)2 - ' ' ' ' 4d4sin(/3:r/2)2 

what is always the case. Summarizing, we obtain that this scheme is unconditionally stable, if 



> -i and unstable otherwise. 



5.2. Integral relations and 2nd order central approximations. Using direct 2nd order 
central approximations for both t and x, we obtain the following scheme: 

= 





{u^t?~Hti+^+i)~^ 


Aft2 ^^J+2 


- '^u-xi + 


^r' 


We denote d := A——, then the polynomial, describing the scheme is 
Aft, 




P = 


= d^T^Tt - T^Tf + (-2^2 + 2)T^Tt 


-7^2; + d^Tt = 


= T^ATt - 


If- 



d'{T, - lyTt, 

with the second expression being the convenient semi-factorized form 

As usual, we use the stability morphism and simplifications. After them, the stability polynomial 
reads then as 52 + (_2+4d2sin^(a/2)).g+l = 0. Denoted := -1+2^2 sin^ (a/2), i.e. 52+26.9+1 = 
and the roots are h ± y/h^ ~ 1. If ^2 > 1^ then one of the roots has modulus bigger, than one. If 
&2 = 1^ the roots are ±1. If ^2 < 1^ the absolute value of both roots equals &2 4- 1 _ 52 _ -y. Hence, 

&2 < 1^ what is satisfied if and only if d < 1 that is — — < — . The same condition is produced 

Aft A 

with the help of CAD in Mathematica. 

Ai 
This scheme is conditionally stable with the condition for the Courant number d — A—— < 1. 

Aft 

5.3. Explicit integration for t and trapezoid rule for x. Using the explicit integration (that 
is, a backward difference) for t and trapezoid rule for x, we get the following scheme. 

^ {u^xi - Hti + ""+2 + 2(« "^f - 'i.u-xi + ^^ "+1) + «r' - 2"r' + - 



4Ai2 



A 



2 



«+2-2<+2 + <+2) = 



The difference scheme polynomial is 

4A2At2 
TIt"^ - 2TlTt + 2T^Tt^ + T^ - 4T^Tt + T^ + 2T^ - 2Ti + 1 tt^{TIT^ - 2T^T^ + T^) 

4A2At2 
Denote d'^ := — . After performing substitutions, we obtain 172 „ 2bg + b = 0, where b = 

(1 + d"^ tan2(a))~^. Its solutions are straightforward: g = b ± \/P — b. If 52 — 6 > 0, we have 6 > 1 
and hence one root is too big. If 52 — 5 < 0, the absolute value of a root is just b'^ + b — b'^ = b, 
what is not bigger than 1. And of course, 6 < 1 is satisfied, since 6=1 + d"^ ta,n^{a))^^ . Hence, 
this scheme is unconditionally stable. With the help of CAD and Mathematica, we arrive to the 
same conclusion. 

5.4. Higher dimensional A-wave equation. One of the crucial advantages of our approach 
and its implementation is the scalability. That is, we employ the algorithms in the very general 
setting. They can be easily modified for the case of more functions (like u) involved. In particular, 
we are able to generate schemes and test them for stability in the higher-dimensional setting. 
Consider the approach from Subsection 15.21 which led us to the conditionally stable scheme. In 
what follows, we apply the same approximations to all the spatial variables. 

Two spatial dimensions. We have Utt — )^^{uxx + Uyy) = 0. The scheme, which we obtain is 
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" ~ ^^2 y'^] + l,k+l ^"j + l.fe+l ^^ "j + l,fe+i; 



X^ X 

/V /^^1 „^)^1 ^^1, /\ 



2 



n+1 Oo,"+1 _L ,,"+1 "l /'o,"+l Oo,"+1 ^ „,"+! 



\"-j+2.k+l ^"i + l.fe+l ^^ "i.fe+l^ A„,2 v"7 + l,fc+2 ^" 



'U„- 



for the Courant numbers dx '■— )^— — , dy :— A- 



/^rj.2 V"'j+2,fe+l ^"-'j + Lk+l ^ "■j.k+ll ^2 V"'j + l,fc+2 ^"j + l.fe+l ^ "^J + l.k) 

In a semi-factorized form, the scheme looks as foUows 

TxTy{Tt - If - dl ■ (T, - ifTyTt - dl ■ Tx{Ty - ifTt = 0. 

The stabihty polynomial in a simplified form is 

g^ - 2(4 cos(^^) + 4 cos{(3y) - 4 - ^ - 2) • .g + 1 = 0. 

Using CAD, we conclude, that this scheme is conditionally stable with the condition 4 + 4 — 1 

At , ._ . Ai 

Ax' ^ Ay 

Three spatial dimensions. The corresponding equation is Uu — )^'^iuxx + Uyy + Uzz) — 0. 
The difference scheme is analogous to the two-dimensional one, in a semi-factorized form it has the 
following form (from which one easily deduces, how the scheme looks in yet higher dimensions): 

TxTyTz{Tt - If - dl ■ {Tx - ifTyTzTt - 4 • Tx{Ty - ifT.Tt ~ dl ■ TxTy{T, - ifTt - 0. 

Running CAD, we obtain, that this scheme, as its lower-dimensional analogues, is conditionally 

stable with the condition dl + dy + dl < 1 for the Courant numbers dx, dy and d^ := ^^7— • 

6. Dispersion Analysis 

6.1. Continuous Dispersion. Recall, that a Fourier node in n-l-1 dimensions is the function 
of the form 



e'ii'^'^)-'^*\ {k,x) := Vfc,x 






Respectively, in 1+1 dimensions it is just Q^ikx-ujt) ^ Qj-^g obtains continuous dispersion from the 
given linear PDE by substituting Fourier nodes into the equation and deriving a relation w = w(fc) 
from the result. 

Example 6.1. For the equation Uu — ^^Uxx — we have 

= ( A _ a2 A)e«('=--"*) = _e^(fe--*) . (^2 _ ;^2^2) 

Hence, uj — ±Xk is the continuous dispersion relation for the X-wave equation. 
We can write down the action of partial derivatives on a Fourier mode. Namely, 



dt- dx 



Hence, the monomial in partial differentiations has its eigenvalue 

^ rr ^(^^^i{k,x)-c.t)■^ ^ ^_^^^a TT(ifcj.)''7 . ^^^i{k,x)^u.t)■^ 

Let us denote F == (,^{k.x)^iujt) _ rj.^^^ d"{F) ^ c{a) ■ F, where a := (a, &i, . . . , &„) G W'+\ 
Extending this action by linearity to the ring of partial differentiations with constant coefficients 
R = K[dt, dxi , . . . , dx„], we are able to compute the eigenvalue of any polynomial from R on F: 

P{F) = Y.p^d'^iF) = (Ep«c(a)) • F. 
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Then, the continuous dispersion relation is obtained by solving with respect to uo the equation 

'^PaC{a) ^0, Pa e K, Ca <E /^(fci , . . . , fc„, w), 

which is called the continuous dispersion equation (CDE) for P. 

Example 6.2. For 1 + n-dimensional heat equation Ut — a^ ■ ^^^i^x^xj = the continuous 
dispersion equation is 



— —iu) — a^ 2_, *^^? "^^^ ^ = —ia^ >, ^? 






Example 6.3. For 1 + n-dimensional modified Xi-wave equation Utt — X]?=i '^? ' "2;ja;i — ^he 



continuous dispersion equation is w ~ iA/X]?=i ^^i^?- 

6.2. Discrete Dispersion. In the discrete case, we consider a discrete Fourier node, correspond- 
ing to the grid point (t,„, {xi)i^,. . . , (x„);„), 

n 

j=i 

One substitutes a discrete Fourier node into the difference scheme and derives a relation uj ~ Lu{k) 
from the result. Let us write down the formula for the eigenvalue of a monomial: 

n n 

As in the continuous case, we extend this action by linearity to polynomials. For a polynomial 
PGif[rt,T,,,...,r,Jonehas 

p{Fn = Y.Pc.T'^{Fn = {Y.v»c{a)) ■ Fr, 

a a 

SO we solve the discrete dispersion equation (DDE) for P, 

^PaC{a) ^0, Pa e K, Ca G K{{kj},Lu), 
a 

and obtain the discrete dispersion relation. Note, that in contrast to the continuous case, this 

relation is not of polynomial form in general. 

Presenting discrete Fourier nodes via trigonometrical functions, we are able to compute discrete 

dispersion relations symbolically. We prefer not to use the de Moivre's formula, but to express 

dispersion relations in terms of sine and cosine of a single argument. 

We work in the commutative ring C{At, Ax)[sint,cost,{sinj,cosj}] modulo the ideal {{sin3 + 

cos^ — 1}, sin^ + coSf ~ 1), where cosj := cos(fcAxj), cost -^ cos(a;Ai). Then, 

n n 

Tt n Tx^i^n - (cos* -*sin,)'^ Y[{cos, +»sin,f ^ • (P/"). 

Example 6.4. Consider the X-wave equation Uu — X'^u^x = and the difference scheme 
d^T^Tt - TxT^ + {-2(f + 2)TxTt - T, + dPTt = 0, 

obtained with the 2"'^ order central approximations for x and t, we denote d := A-^— . 

Performing computations, we obtain after simplification d^ cos^ —cost -\- 1 — d^ = 0, that is 
cos(cjAt) = 1 — d^{l — cos(fcAx)). In the stability limit d — > 1, w;e have cos{ujAt) = cos{kAx), 
hence w = ±-^fc + 2T:m,m G Z. Since d — > 1 implies -^ — > A, in the stability limit the discrete 
dispersion relation becomes cj — ±Xk + 2TTm, where for m ^ we recover the continuous dispersion 
relation. 
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7. Conclusion and Future Work 



The advantages of the methods we propose include, among other, their scalabihty and tendency 
towards automatization. Indeed, we do not make distinction between classical types of PDE's (hy- 
perbolic, elliptic, parabolic). Thus these methods are very general. Symbolic methods are able to 
generate automatically many difference schemes of standard linear PDE's, as it was demonstrated 
in [11] and by ourselves. 

An important issue for the future research is the (partial) algebraization of the consistency of 
a generated scheme with the differential equation. Provided such a check, one could work with 
general multi-parametric schemes, where the conditions on parameters arise from the consistency 
check and the symbolic stability approach. 

We decided not to include the treatment of systems of linear PDE's in this paper. However, we 
want to remark, that by the rewriting system approach the number of the discretized equations is 
exactly the number of PDE's one started with. By using Grobner bases, we get in general more 
equations, which, in turn, reveal the interplay between discretized equations. Such interplay is 
not detected by the rewriting approach at all; it seems to us that numerists just ignored such 
interplay. Hence, this issue need to be investigated further. 

Christian Dingier (TU Kaiserslautern, Germany) has been developing a new package for Singular 
with QEPCAD as an engine for cylindrical algebraic decomposition. This package extends the 
tools for the generation of finite difference schemes to the cases of a single linear PDE and of a 
system of linear PDE's. Another problem for the further research is the generalization of von 
Neumann stability for systems, which is clear only for some classes of equations. 
The generalization of the methods we described goes in several different directions: allowing 
variables coefficients and/or allowing nonlinearity of expressions in functions. The generation of 
schemes in such settings is still possible (by the rewriting system approach), but even the notion 
of stability, to the best of our knowledge, needs to be investigated depending on the particular 
class of equations. 

Very important question is the role of differential and difference Grobner bases for nonlinear 
equations in the scheme generation and analysis. The recent works [111 110] show, that in some 
cases the systematic use of interplay between equations can produce more unviversal (though, of 
course, more complicated) schemes, which hints at big potential of these methods. 
In connection with this, a new theory for infinite (but still constructively approachable!) difference 
Grobner bases, arising from the letterplace philosophy [17], is of great interest. 
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Appendix. Tools for finite difference schemes in Singular 

8. An introduciton to the system Singular and its language 

Next we describe shortly by examples how to read Singular language and how to obtain and 
interpret the output - as far as it is used to generate a difference scheme. More details are found 
in the Singular manual at www.singular.uni-kl.de. 

• Defining an algebra. At first a domain for the computation has to be defined. This is done 
by the following input: 

ring R = (0,ro ,K,dt ,dh) , (Tx,Tt) , (c,dp) ; 
It defines a polynomial ring, denoted by R, which represents an algebra over a prime field. 
The first brackets defines the coefficient field. It starts with a number for the character- 
istic, here stands for the prime field of rational numbers, which is always used in our 
considerations. A list of names for parameters follows. These parameters are constants 
- as step size in space or time {dh,dt) or constants of the equations {ro,K). Hence we 
have defined K as field of coefficients, here: rational functions over Q in the parameter 
ro, K, dt, dh. The list of the second brackets corresponds to names of the variables Tx, Tt, 
which stand for commuting polynomial variables over the coefficient field K (we restrict us 
for this description to the commutative case only). The names are almost arbitrary strings, 
hence names can be used to indicate the meaning - here the shift operators. The string 
in the third brackets explains the monomial order used for Grobner basis computations: 
dp - degree reverse lexicographical ordering for polynomials. A small c at first place sort 
polynomial vectors by components first in descending order, i.e., ei > 62 > . . ., then by 
the monomial ordering. 

• Evaluating the constants: 

Wanting to evaluate the constants in a resulting expression you have to create a new ring, 
a mapping and substitutions: 

ring r = (0,dt ,dh) , (Tx,Tt ,ro,K) , (c,dp) ; 

def Dbl = imap(R,Db) ; 

number num_ro = dh/2; number num_K = 4*dt~2 - dh~2; // constants 

def Db2 = subst (Dbl ,ro,num_ro) ; 
Db2 = subst ( Db2, K, number_K) ; 
Here 'Ob' is the name of the type (polynomial, matrix or ideal) you are considering, 'def 
Obi ' creates an object of type of the right hand side. A new ring is introduced because a 
substitution is defined only for variables. 

• Creating a matrix: Starting with a linear system of PDE's, adding approximation rules, we 
end up with an extended system AU — 0, see above. We need only the matrix with entries 
in the ring R of shift operators: 

matrix A [3] [3] = 

(-Tx*Tt~2+Tx) , (Tx~2*Tt - Tt) , , 

0, (dh/2)*(Tx+l), 1-Tx, 

(dt/2)*(Tt+l), 0, 1-Tt; 
In the definition of a matrix you have to indicate row- and column-size, and on the right 
hand side just the list of polynomials. 

• Elimination of components: 

We have to eliminate all components of the vector U that stand for a differential operator, 
corresponding to Gauss-elimination over R with the matrix. In this example the anony- 
mous vector U stands for {uitUxtuY ■ We want to produce with A a row, having entries 
only in the last component (s). This is done most efficiently by a Grobner basis compu- 
tation of the module generated by the columns of the transposed matrix in the indicated 
monomial ordering. The last non-zero component(s) of the first generator(s) correspond to 
the difference scheme, here: 

module M = transpose (A) ; 

module Ml = std(M) ; 



22 V. LEVANDOVSKYY AND B. MARTIN 

print (Ml) ; 

Ml [3,1]; 
The command 'print ' returns only the first digits of the string corresponding to any entry. 
Type the indices of an entry in square brackets, to get it in fuh length. Here the resulting 
equation is returned as: 

(-dt)*Tx~2*Tt+(dh)*Tx*Tt~2+(2*dt-2*dh)*Tx*Tt+(dh)*Tx+(-dt)*Tt 

8.1. Tools for diflference schemes. The library findifs.lib has been created to automate 
numerous processes, taking place while working with finite difference schemes. An important role 
is played by sophisticated routines, transforming the different forms of objects into some classical 
ones. One can generate complicated schemes and easily present them in, say, nodal form. At the 
same time, one keeps the polynomial operator presentation, which is used in e.g. stability analysis. 

• decoef (P,n) ; P a poly, n a number. 

Decomposes the poly P into summands with respect to the presence of a number n in the 
coefficients, returns an ideal in usually two generators. 
Example: 

ring r = (0,dh,dt) , (Tx,Tt) ,dp; 

poly P = (4*dh"2-dt)*Tx'3*Tt + dt*dh*Tt~2 + dh*Tt ; 

P; 

==> (4*dh~2-dt)*Tx~3*Tt+(dh*dt)*Tt~2+(dh)*Tt 

decoef (P,dt) ; 
==> _[l] = (4*dh~2)*Tx"3*Tt+(dli)*Tt // the part, not containing dt 
_[2]=(-dt)*Tx~3*Tt+(dh*dt)*Tt~2 // the part which contains dt 
decoef (P,dh) ; 
==> _ [1] =(-dt)*Tx~3*Tt // the part, not containing dh 

_ [2] = (4*dh~2) *Tx"3*Tt+ (dh*dt ) *Tt '2+ (dh) *Tt 

• difpoly2texS,P [,Q] ); S an ideal, P and optional Q are lists. 

Presents the difference scheme, given in the ideal S, in the nodal form. The ideal S thought 
to be the result of decoef, list P contains parameters, which will be controlled in order to 
remain in numerators. The optional list Q contains polynomials, which will be added to 
the scheme (written in the function u) the part in terms of a function p. 
Example: 

ring r = (0,dh,dt ,V) , (Tx,Tt) ,dp; 
poly M = (2*dh*Tx+dt)'2*(Tt-l) + V*Tt*Tx; 
M; 
==> (4*dh~2)*Tx~2*Tt+(-4*dh~2)*Tx~2+(4*dh*dt+V)*Tx*Tt+(-4*dh*dt)*Tx+(dt~2)*Tt+(-dt~2) 
ideal 1 = decoef (M,dt) ; 

I; 

==> 1 [1] = (4*dh-2) *Tx-2*Tt+ (-4*dh-2) *Tx-2+ (V) *Tx*Tt 

l[2]=(4*dh*dt)*Tx*Tt+(-4*dh*dt)*Tx+(dt~2)*Tt+(-dt~2) 
list L; L[l] = V; 
difpoly2tex(l ,L) ; 
==> \frac{l}{4 \tri t}\cdot (u-{n+l}_{j+2}-u-{n}_{j+2}+\f rac{ \nu}{4 \tri h "{2}} 
u-{n+l}_{j+l})+ \frac{l}{4 \tri h}\cdot (u-{n+l}_{j+l}-u-{n}_{j+l}+ 
\frac{ \tri t}{4 \tri h} u~{n+l}_{j}+\f rac{- \tri t}{4 \tri h} u~{n}_{j» 

The last output, compiled with TcX, produces 

4At ^"^+2 "j+2+4^/^2"j+ii + 4A/, ^"^+1 "^■+i + 4A/i ^' ^AAh^^'' 
Now let us illustrate the use of the optional list Q. 

ring D = (0,ro,K,dt ,dh) , (Tx,Tt) , (c,Dp) ; 
poly U = (-K*dt)*Tx~2*Tt+(K*dt)*Tt; 
poly P = (-2*ro*dh)*Tx*Tt+(2*ro*dh)*Tx; 
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list V; V[l] = K; V [2] = ro ; 

difpoly2tex(-U,V,-P) ; // we want P to be in terms of p, and U in terms of u 
==> \frac{K}{2 \tri h}\cdot (u~{n+l}_{j+2}-u"{n+l}_{j})+ \frac{ \rho}{ \tri t} 
\cdot (p-{n+l}_{j+l}-p-{n}_{j+l» 
That is, we have produced the nodal form of scheme for two functions u and p: 
K n 

2Ah ^ ■''+^ ^' ' At ^^+^ P]+i)- 
exp2pt (P [,L] ) ; convert a polynomial M into the TeX format, in nodal form 
mon2pt (P [,L] ) ; convert a monomial M into the TeX format, in nodal form 
texcoef (n) ; converts the number n into TeX 

npar(n) ; search for 'n' among the parameters and returns its number 
replace (s, what , with) ; replaces in s all the substrings with a given string 
xchange(w,a,b) ; exchanges two substrings of a string 
par2tex(s) ; converts special characters to TeX in s. 
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